I.I: Manually build graphs from adjacency matrices.

Figure 1. Cartoon of graph we will build manually via its adjacency matrix.

adj_mat = rbind(c(0,3,1,0,0,0,0,0,0), #A's neighbors
                c(3,0,2,2,0,0,0,0,0), #B's neighbors
                c(1,2,0,0,2,0,0,0,0), #C's neighbors
                c(0,2,0,0,1,0,1,1,0), #D's neighbors
                c(0,0,2,1,0,2,0,2,0), #E's neighbors
                c(0,0,0,0,2,0,0,0,0), #F's neighbors
                c(0,0,0,1,0,0,0,1,0), #G's neighbors
                c(0,0,0,1,2,0,1,0,1), #H's neighbors
                c(0,0,0,0,0,0,0,1,0) #I's neighbors
                )
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)

## IGRAPH 419525a UNW- 9 12 -- 
## + attr: name (v/c), weight (e/n)
## + edges from 419525a (vertex names):
##  [1] A--B A--C B--C B--D C--E D--E D--G D--H E--F E--H G--H H--I

I.II: Learn a graph from data.

# Load the Miller2015_Heparin dataset
data(Miller2015_Heparin)
data(heparin_metabolites)
data = Miller2015_Heparin
# Only include metabolites that are present in >90% reference samples.
fil.rate = anno$`Times identifed in all 200 samples`/200
names(fil.rate) = rownames(anno)
data = data[which(fil.rate>0.90), ]
#dim(data)
# Remove any metabolites where any profile has a z-score > 1000. These are likely imputed raw values that were not z-scored.
rmMets = names(which(apply(data, 1, function(i) any(i>20))))
if (length(rmMets)>0) {
  data = data[-which(rownames(data) %in% rmMets),]
}
#dim(data)

# Get the diagnoses associated with the samples
data(Miller2015_diags)
# Get data from all patients with Argininemia
arg_data = data[,which(diagnoses$diagnosis=="Argininemia")]
# Add surrogate disease and surrogate reference profiles based on 1 standard deviation around profiles from real patients to improve rank of matrix when learning Gaussian Markov Random Field network on data.
arg_data = data.surrogateProfiles(arg_data, 1, FALSE, TRUE, ref_data = data[,which(diagnoses$diagnosis=="No biochemical genetic diagnosis")])
#dim(arg_data)

# Learn a Gaussian Markov Random Field model using the Graphical LASSO in the R package "huge". 
# Select the regularization parameter based on the "STARS" stability estimate.
# require(huge)
#This will take 30 seconds - 1 minute.
arg = huge(t(arg_data), method="glasso")
plot(arg)
# This will take 3-5-ish minutes.
arg.select = huge.select(arg, criterion="stars")
plot(arg.select)
# This is the regularization parameter the STARS method selected.
print(arg.select$opt.lambda)
# This is the corresponding inverse of the covariance matrix that corresponds to the selected regularization level.
arg_icov = as.matrix(arg.select$opt.icov)
# Remove all "self" edges, as we are not interested in self-relationships.
diag(arg_icov) = 0
rownames(arg_icov) = rownames(arg_data)
colnames(arg_icov) = rownames(arg_data)
# Convert adjacency matrices to an igraph object.
ig_arg = graph.adjacency(arg_icov, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig_arg)

II.I From a starting node.

Run the following code, then go to the directory, and open all diffusionEventMovie*.png files all at once. Starting from the first image, view how the probability diffusion algorithm works to diffuse 100% probability to the rest of the graph. Be sure to pay attention to the recursion level listed in the title of each image, to imagine where in the call stack the algorithm is at the captured time the image was generated.

# Set some global parameters for the Probability Diffusion Algorithm. 
p0=0.0
p1=1.0
thresholdDiff=0.01

G=vector(mode="list", length=length(V(ig)$name))
G[1:length(G)] = 0
names(G) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
startNode = "A"
visitedNodes = startNode
# Probability diffusion truncates at
thresholdDiff=0.01
coords = layout.fruchterman.reingold(ig)
V(ig)$x = coords[,1]
V(ig)$y = coords[,2]
# Global variable imgNum
imgNum=1
G_new = graph.diffuseP1Movie(p1, startNode, G, visitedNodes, ig, recursion_level=1, output_dir = "/Users/lillian.rosa/Downloads/CTD/vignette/diffusion_event")
# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))

## [1] 1

# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]

## $B
## [1] 0.40625

II.II From a starting node, after visiting previous nodes.

Now, delete all diffusionEventMovie*.png files from your current directory, and run the following code. View the new image stack in the same way we did previously.

# Now let's see how the probability diffusion algorithm diffuses probability after B has been "stepped" into.
visitedNodes = c("A", "B")
startNode = "B"
imgNum=1

G_new = graph.diffuseP1Movie(p1, startNode, G, visitedNodes, ig, 1, output_dir = "/Users/lillian.rosa/Downloads/CTD/vignette/diffusion_event2")
# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))

## [1] 1

# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]

## $D
## [1] 0.2791667

II.III Diffusing through visited nodes, based on connectivity.

Sometimes the startNode is “stranded” by a bunch of visited nodes. The diffusion algorithm diffuses “through” visited nodes, so that nodes in the same connected component can be prioritized over nodes in a different connected component, or “island nodes” (e.g. in the below code snippet, “I” is an island node). This only works currently for nodes 2 hops away from the current startNode, however.

adj_mat = rbind(c(0,1,2,0,0,0,0,0,0), #A's neighbors
                c(1,0,3,0,0,0,0,0,0), #B's neighbors
                c(2,3,0,0,1,0,0,0,0), #C's neighbors
                c(0,0,0,0,0,0,1,1,0), #D's neighbors
                c(0,0,1,0,0,1,0,0,0), #E's neighbors
                c(0,0,0,0,1,0,0,0,0), #F's neighbors
                c(0,0,0,1,0,0,0,1,0), #G's neighbors
                c(0,0,0,1,0,0,1,0,0), #H's neighbors
                c(0,0,0,0,0,0,0,0,0) #I's neighbors
                )
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)

## IGRAPH f620114 UNW- 9 8 -- 
## + attr: name (v/c), weight (e/n)
## + edges from f620114 (vertex names):
## [1] A--B A--C B--C C--E D--G D--H E--F G--H

adjacency_matrix = list(adj_mat)

# Now let's see how the probability diffusion algorithm diffuses probability after B has been "stepped" into "C" and then "A". As you can see, startNode "A" is surrounded by visited nodes "B" and "C". It needs to be smart enough to weigh "E" and "F" before "D", "H", "G" and "I".
visitedNodes = c("B", "C", "A")
startNode = "A"
G_new = graph.diffuseP1(1.0, startNode, G, visitedNodes, 1, verbose=TRUE)

## [1] "      prob. to diffuse:1.000000 startNode: A, visitedNodes: B, C, A"
## [1] "      Weighted_edges.sum=1.000000"
## [1] "      added diffused  probability 1.000000 to child #1: E"
## [1] "      subtracted 0.500000 from startNode's neighbor #1: E and sent to its own neighbors."

# Inherited probabilities across all nodes should add to 1.
sum(unlist(G_new))

## [1] 1

# Which node inherited the highest probability from startNode?
G_new[which.max(G_new)]

## $E
## [1] 0.5

III.I Adaptive (“With Memory”) Network Walker

adj_mat = rbind(c(0,3,1,0,0,0,0,0,0), #A's neighbors
                c(3,0,2,2,0,0,0,0,0), #B's neighbors
                c(1,2,0,0,2,0,0,0,0), #C's neighbors
                c(0,2,0,0,1,0,1,1,0), #D's neighbors
                c(0,0,2,1,0,2,0,2,0), #E's neighbors
                c(0,0,0,0,2,0,0,0,0), #F's neighbors
                c(0,0,0,1,0,0,0,1,0), #G's neighbors
                c(0,0,0,1,2,0,1,0,1), #H's neighbors
                c(0,0,0,0,0,0,0,1,0) #I's neighbors
                )
rownames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
colnames(adj_mat) = c("A", "B", "C", "D", "E", "F", "G", "H", "I")
# Convert adjacency matrices to igrpah objects for all three graphs.
ig = graph.adjacency(adj_mat, mode="undirected", weighted=TRUE, add.colnames = "name")
print(ig)

## IGRAPH 5ee3d84 UNW- 9 12 -- 
## + attr: name (v/c), weight (e/n)
## + edges from 5ee3d84 (vertex names):
##  [1] A--B A--C B--C B--D C--E D--E D--G D--H E--F E--H G--H H--I

adjacency_matrix = list(adj_mat)

perms = mle.getPermMovie_memory(subset.nodes = c("A", "B"), ig, output_filepath = "/Users/lillian.rosa/Downloads/CTD/vignette/biased_walker", movie=TRUE, zoomIn = FALSE)

## [1] "Calculating permutation 1 of 2..."

## [1] "Calculating permutation 2 of 2..."

# Get perms as list object, with no images generated
perms = mle.getPerms_memory(S = c("A", "B"), G)

## [1] "Calculating permutation 1 of 2."
## [1] "Calculating permutation 2 of 2."

III.II Unbiased (“Memoryless”) Network Walker

# This network walker is less sensitive, but node permutations can be calculated ahead of time (along the lines of a dynamic programming solution), 
#      so thousands of subsets can be reasoned about very quickly.

# Generate PNGs to animate the network walker.
perms = mle.getPermMovie_memoryless(subset.nodes = c("A", "B"), ig, output_filepath = "/Users/lillian.rosa/Downloads/CTD/vignette/unbiased_walker", movie=TRUE, zoomIn = FALSE)

## [1] "Calculating permutation 1 of 2..."

## [1] "Calculating permutation 2 of 2..."

# Get perms as list object, with no images generated
S = c("A", "B")

perms = list()
for (n in 1:length(S)) {
  perms[[n]] = mle.getPermN_memoryless(n, G, S, misses.thresh = log2(length(G)))
}

## [1] "Calculating permutation 1 of 9."
## [1] "Calculating permutation 2 of 9."

names(perms) = S

IV.II Get k node permutations.

# Get the adaptive network walker's paths starting from each node in the subset S_arg.
perms = mle.getPerms_memory(S = names(S_arg), G)

## [1] "Calculating permutation 1 of 15."
## [1] "Calculating permutation 2 of 15."
## [1] "Calculating permutation 3 of 15."
## [1] "Calculating permutation 4 of 15."
## [1] "Calculating permutation 5 of 15."
## [1] "Calculating permutation 6 of 15."
## [1] "Calculating permutation 7 of 15."
## [1] "Calculating permutation 8 of 15."
## [1] "Calculating permutation 9 of 15."
## [1] "Calculating permutation 10 of 15."
## [1] "Calculating permutation 11 of 15."
## [1] "Calculating permutation 12 of 15."
## [1] "Calculating permutation 13 of 15."
## [1] "Calculating permutation 14 of 15."
## [1] "Calculating permutation 15 of 15."

IV.III Convert to bitstrings.

ptBSbyK = mle.getPtBSbyK_memory(names(S_arg), perms)

IV.IV Get encoding length of minimum length codeword.

ptID = colnames(arg_data)[1]
data_mx.pvals=apply(arg_data[,which(colnames(arg_data) %in% diagnoses$id)], c(1,2), function(i) 2*pnorm(abs(i), lower.tail=FALSE))
res = mle.getEncodingLength(ptBSbyK, t(data_mx.pvals), ptID, G)

IV.V Get probability of node subset.

2^-res[,"IS.alt"]

##  [1] 1.615509e-03 8.077544e-04 4.038772e-04 2.019386e-04 1.009693e-04
##  [6] 5.048465e-05 2.524233e-05 7.888227e-07 7.703347e-10 1.925837e-10
## [11] 9.629183e-11 4.814592e-11 2.407296e-11 1.203648e-11 6.018240e-12

IV.V Get p-value of variable length encoding vs. fixed length encoding.

2^-res[,"d.score"]

##  [1] 1.000000e+00 6.474570e-03 6.292531e-05 8.176575e-07 1.329077e-08
##  [6] 2.597783e-10 5.934544e-12 2.481370e-12 3.743036e-11 2.454003e-12
## [11] 8.864218e-14 3.498948e-15 1.498843e-16 6.929281e-18 3.435771e-19